# Replication code for "Wartime Violence and Post-War Women's Representation -- Online Appendix (OA)"

# clear workspace
rm(list=ls())

# install relevant packages before loading them 
# install.packages("multiwayvcov")
# install.packages("lmtest")
# install.packages("coin")
# install.packages("betareg")
# install.packages("car")
# install.packages("censReg")
# install.packages("MASS")

# load packages
library(multiwayvcov)
library(lmtest)
library(coin)
library(betareg)
library(car)
library(censReg)
library(MASS)

# ensure that you have set the appropriate working directory
# setwd("...")

dataCandidates <- read.csv("Candidates_Data.csv")
dataCouncilors <- read.csv("Councilors_Data.csv")
dataCandidates <- dataCandidates[1:87,]

##### OA2: Descriptive Statistics #####
# Table OA2.1: Descriptive Statistics
# Candidates
round(mean(dataCandidates$Candidates),2) # Mean
round(median(dataCandidates$Candidates),2) # Median
round(sd(dataCandidates$Candidates),2) # St. dev.
round(min(dataCandidates$Candidates),2) # Min.
round(max(dataCandidates$Candidates),2) # Max.

# Councilors
#   for 1990 election (1990 Councilors)
round(mean(dataCandidates$Y1990_councilors),2) # Mean
round(median(dataCandidates$Y1990_councilors),2) # Median
round(sd(dataCandidates$Y1990_councilors),2) # St. dev.
round(min(dataCandidates$Y1990_councilors),2) # Min.
round(max(dataCandidates$Y1990_councilors),2) # Max.

#   for 2000 election
round(mean(dataCandidates$Councilors_2000),2) # Mean
round(median(dataCandidates$Councilors_2000),2) # Median
round(sd(dataCandidates$Councilors_2000),2) # St. dev.
round(min(dataCandidates$Councilors_2000),2) # Min.
round(max(dataCandidates$Councilors_2000),2) # Max.

#   for 1990 and 2000 elections
Total_Councilors <- c(dataCandidates$Y1990_councilors, dataCandidates$Councilors_2000)
round(mean(Total_Councilors),2) # Mean
round(median(Total_Councilors),2) # Median
round(sd(Total_Councilors),2) # St. dev.
round(min(Total_Councilors),2) # Min.
round(max(Total_Councilors),2) # Max.

# Casualty
round(mean(dataCandidates$Casualty),2) # Mean
round(median(dataCandidates$Casualty),2) # Median
round(sd(dataCandidates$Casualty),2) # St. dev.
round(min(dataCandidates$Casualty),2) # Min.
round(max(dataCandidates$Casualty),2) # Max.

# log(Casualty)
round(mean(dataCandidates$log_Casualty),2) # Mean
round(median(dataCandidates$log_Casualty),2) # Median
round(sd(dataCandidates$log_Casualty),2) # St. dev.
round(min(dataCandidates$log_Casualty),2) # Min.
round(max(dataCandidates$log_Casualty),2) # Max.

# Confirmed dead
round(mean(dataCandidates$Confirmed_dead),2) # Mean
round(median(dataCandidates$Confirmed_dead),2) # Median
round(sd(dataCandidates$Confirmed_dead),2) # St. dev.
round(min(dataCandidates$Confirmed_dead),2) # Min.
round(max(dataCandidates$Confirmed_dead),2) # Max.

# log(Confirmed dead)
round(mean(dataCandidates$log_Confirmed_dead),2) # Mean
round(median(dataCandidates$log_Confirmed_dead),2) # Median
round(sd(dataCandidates$log_Confirmed_dead),2) # St. dev.
round(min(dataCandidates$log_Confirmed_dead),2) # Min.
round(max(dataCandidates$log_Confirmed_dead),2) # Max.

# Urban share
round(mean(dataCandidates$Urban_share),2) # Mean
round(median(dataCandidates$Urban_share),2) # Median
round(sd(dataCandidates$Urban_share),2) # St. dev.
round(min(dataCandidates$Urban_share),2) # Min.
round(max(dataCandidates$Urban_share),2) # Max.

# Income per capits
round(mean(dataCandidates$Income_per_capita),2) # Mean
round(median(dataCandidates$Income_per_capita),2) # Median
round(sd(dataCandidates$Income_per_capita),2) # St. dev.
round(min(dataCandidates$Income_per_capita),2) # Min.
round(max(dataCandidates$Income_per_capita),2) # Max.

# Ethnic polarization
round(mean(dataCandidates$Ethnic_polarization),2) # Mean
round(median(dataCandidates$Ethnic_polarization),2) # Median
round(sd(dataCandidates$Ethnic_polarization),2) # St. dev.
round(min(dataCandidates$Ethnic_polarization),2) # Min.
round(max(dataCandidates$Ethnic_polarization),2) # Max.

# Share of women
round(mean(dataCandidates$Share_of_women),2) # Mean
round(median(dataCandidates$Share_of_women),2) # Median
round(sd(dataCandidates$Share_of_women),2) # St. dev.
round(min(dataCandidates$Share_of_women),2) # Min.
round(max(dataCandidates$Share_of_women),2) # Max.

# Figure OA2.1: Variable Density Plots
plot(density(dataCandidates$Candidates), main='Candidates')
plot(density(dataCandidates$Y1990_councilors), main='Councilors (in 1990)')
plot(density(dataCandidates$Councilors_2000), main='Councilors (in 2000)')
plot(density(Total_Councilors), main='Councilors (in 1990 and 2000)')
plot(density(dataCandidates$Casualty), main='Casualty')
plot(density(dataCandidates$log_Casualty), main='log(Casualty)')
plot(density(dataCandidates$Confirmed_dead), main='Confirmed dead')
plot(density(dataCandidates$log_Confirmed_dead), main='log(Confirmed dead)')
plot(density(dataCandidates$Urban_share), main='Urban share')
plot(density(dataCandidates$Income_per_capita), main='Income per capita')
plot(density(dataCandidates$Ethnic_polarization), main='Ethnic polarization')
plot(density(dataCandidates$Share_of_women), main='Share of women')

#############################################################################
# PLEASE SEE "Heat_Maps.zip" TO REPLICATE "Figure OA2.2: Variable Heat Maps" 
#############################################################################

##### OA3: Missing Data #####
data101 <- read.csv("Candidates_Data.csv")

### NOTE ABOUT RESULTS FROM TABLE OA3.1: please note that models presented in Table OA3.1 test a directional hypothesis (violence has a
### positive association with women's participation) and employ one-tailed t-tests. Therefore, in order to obtain the correct
### p-values, (a) divide the p-values presented in the output by 2 and (b) remember that only positive coefficient estimates
### can be marked as statistically significant.

# Table OA3.1: Wartime Violence and Female Candidates, Not Controlling for 1990 Councilors
# Model 1: Candidates
summary(mod1 <- lm(Candidates ~ log_Casualty + Urban_share + Income_per_capita + Ethnic_polarization +
                     Share_of_women, data=data101))
nobs(mod1)

# Model 2: Candidates
summary(mod2 <- lm(Candidates ~ log_Confirmed_dead + Urban_share + Income_per_capita + Ethnic_polarization +
                     Share_of_women, data=data101))
nobs(mod2)

dataMissing <- read.csv("Missing_Muns.csv")
round(mean(dataCandidates$Urban_share),0) # 34 average urbanization for present municipalities
round(mean(dataMissing$Urban_share),0) # 27 average urbanization for exlcluded municipalities

# Table OA3.2: Wilcoxon Tests Comparing Municipalities in the Sample to Municipalities with Missing Data
# Candidates
# # Missing
sum(!complete.cases(dataMissing$Candidates))
# # in Comparison Group
sum(complete.cases(dataMissing$Candidates))
# p-value
wilcox.test(dataCandidates$Candidates, dataMissing$Candidates) 

# Councilors (in 1990)
# # Missing
sum(!complete.cases(dataMissing$Y1990_councilors))
# # in Comparison Group
sum(complete.cases(dataMissing$Y1990_councilors))
# p-value
wilcox.test(dataCandidates$Y1990_councilors, dataMissing$Y1990_councilors)

# Councilors (in 2000)
# # Missing
sum(!complete.cases(dataMissing$Councilors_2000))
# # in Comparison Group
sum(complete.cases(dataMissing$Councilors_2000))
# p-value
wilcox.test(dataCandidates$Councilors_2000, dataMissing$Councilors_2000)

# Casualty
# # Missing
sum(!complete.cases(dataMissing$Casualty))
# # in Comparison Group
sum(complete.cases(dataMissing$Casualty))
# p-value
wilcox.test(dataCandidates$Casualty, dataMissing$Casualty)

# log(Casualty)
# # Missing
sum(!complete.cases(dataMissing$log_Casualty))
# # in Comparison Group
sum(complete.cases(dataMissing$log_Casualty))
# p-value
wilcox.test(dataCandidates$log_Casualty, dataMissing$log_Casualty)

# Confirmed dead
# # Missing
sum(!complete.cases(dataMissing$Confirmed_dead))
# # in Comparison Group
sum(complete.cases(dataMissing$Confirmed_dead))
# p-value
wilcox.test(dataCandidates$Confirmed_dead, dataMissing$Confirmed_dead)

# log(Confirmed dead)
# # Missing
sum(!complete.cases(dataMissing$log_Confirmed_dead))
# # in Comparison Group
sum(complete.cases(dataMissing$log_Confirmed_dead))
# p-value
wilcox.test(dataCandidates$log_Confirmed_dead, dataMissing$log_Confirmed_dead)

# Urban share
# # Missing
sum(!complete.cases(dataMissing$Urban_share))
# # in Comparison Group
sum(complete.cases(dataMissing$Urban_share))
# p-value
wilcox.test(dataCandidates$Urban_share, dataMissing$Urban_share)

# Income per capita
# # Missing
sum(!complete.cases(dataMissing$Income_per_capita))
# # in Comparison Group
sum(complete.cases(dataMissing$Income_per_capita))
# p-value
wilcox.test(dataCandidates$Income_per_capita, dataMissing$Income_per_capita)

# Ethnic polarization
# # Missing
sum(!complete.cases(dataMissing$Ethnic_polarization))
# # in Comparison Group
sum(complete.cases(dataMissing$Ethnic_polarization))
# p-value
wilcox.test(dataCandidates$Ethnic_polarization, dataMissing$Ethnic_polarization)

# Share of women
# # Missing
sum(!complete.cases(dataMissing$Share_of_women))
# # in Comparison Group
sum(complete.cases(dataMissing$Share_of_women))
# p-value
wilcox.test(dataCandidates$Share_of_women, dataMissing$Share_of_women)

data101Candidates <- c(dataCandidates$Candidates, dataMissing$Candidates) 
data101Y1990_councilors <- c(dataCandidates$Y1990_councilors, dataMissing$Y1990_councilors) 
data101Councilors_2000 <- c(dataCandidates$Councilors_2000, dataMissing$Councilors_2000) 
data101Casualty <- c(dataCandidates$Casualty, dataMissing$Casualty) 
data101log_Casualty <- c(dataCandidates$log_Casualty, dataMissing$log_Casualty) 
data101Confirmed_dead <- c(dataCandidates$Confirmed_dead, dataMissing$Confirmed_dead) 
data101log_Confirmed_dead <- c(dataCandidates$log_Confirmed_dead, dataMissing$log_Confirmed_dead) 
data101Urban_share <- c(dataCandidates$Urban_share, dataMissing$Urban_share) 
data101Income_per_capita <- c(dataCandidates$Income_per_capita, dataMissing$Income_per_capita) 
data101Ethnic_polarization <- c(dataCandidates$Ethnic_polarization, dataMissing$Ethnic_polarization) 
data101Share_of_women <- c(dataCandidates$Share_of_women, dataMissing$Share_of_women) 
# Table OA3.3: Wilcoxon Tests Comparing Municipalities in the Sample to Filled-Out Data
# Candidates
# # Missing
sum(!complete.cases(data101Candidates))
# # in Comparison Group
sum(complete.cases(data101Candidates))
# p-value
wilcox.test(dataCandidates$Candidates, data101Candidates) 

# Councilors (in 1990)
# # Missing
sum(!complete.cases(data101Y1990_councilors))
# # in Comparison Group
sum(complete.cases(data101Y1990_councilors))
# p-value
wilcox.test(dataCandidates$Y1990_councilors, data101Y1990_councilors)

# Councilors (in 2000)
# # Missing
sum(!complete.cases(data101Councilors_2000))
# # in Comparison Group
sum(complete.cases(data101Councilors_2000))
# p-value
wilcox.test(dataCandidates$Councilors_2000, data101Councilors_2000)

# Casualty
# # Missing
sum(!complete.cases(data101Casualty))
# # in Comparison Group
sum(complete.cases(data101Casualty))
# p-value
wilcox.test(dataCandidates$Casualty, data101Casualty)

# log(Casualty)
# # Missing
sum(!complete.cases(data101log_Casualty))
# # in Comparison Group
sum(complete.cases(data101log_Casualty))
# p-value
wilcox.test(dataCandidates$log_Casualty, data101log_Casualty)

# Confirmed dead
# # Missing
sum(!complete.cases(data101Confirmed_dead))
# # in Comparison Group
sum(complete.cases(data101Confirmed_dead))
# p-value
wilcox.test(dataCandidates$Confirmed_dead, data101Confirmed_dead)

# log(Confirmed dead)
# # Missing
sum(!complete.cases(data101log_Confirmed_dead))
# # in Comparison Group
sum(complete.cases(data101log_Confirmed_dead))
# p-value
wilcox.test(dataCandidates$log_Confirmed_dead, data101log_Confirmed_dead)

# Urban share
# # Missing
sum(!complete.cases(data101Urban_share))
# # in Comparison Group
sum(complete.cases(data101Urban_share))
# p-value
wilcox.test(dataCandidates$Urban_share, data101Urban_share)

# Income per capita
# # Missing
sum(!complete.cases(data101Income_per_capita))
# # in Comparison Group
sum(complete.cases(data101Income_per_capita))
# p-value
wilcox.test(dataCandidates$Income_per_capita, data101Income_per_capita)

# Ethnic polarization
# # Missing
sum(!complete.cases(data101Ethnic_polarization))
# # in Comparison Group
sum(complete.cases(data101Ethnic_polarization))
# p-value
wilcox.test(dataCandidates$Ethnic_polarization, data101Ethnic_polarization)

# Share of women
# # Missing
sum(!complete.cases(data101Share_of_women))
# # in Comparison Group
sum(complete.cases(data101Share_of_women))
# p-value
wilcox.test(dataCandidates$Share_of_women, data101Share_of_women)

##### OA5: Violence and District Magnitude #####
# Table OA5.1: Correlations between Violence and District Magnitude
# Casualty and 1990 Magnitude
cor.test(dataCandidates$Casualty, dataCandidates$Magnitude_1990)
# Casualty and 2000 Magnitude
cor.test(dataCandidates$Casualty, dataCandidates$Magnitude_2000)
# Casualty and Change in Magnitude
cor.test(dataCandidates$Casualty, dataCandidates$Magnitude_Change)

# log(Casualty) and 1990 Magnitude
cor.test(dataCandidates$log_Casualty, dataCandidates$Magnitude_1990)
# log(Casualty) and 2000 Magnitude
cor.test(dataCandidates$log_Casualty, dataCandidates$Magnitude_2000)
# log(Casualty) and Change in Magnitude
cor.test(dataCandidates$log_Casualty, dataCandidates$Magnitude_Change)

# Confirmed dead and 1990 Magnitude
cor.test(dataCandidates$Confirmed_dead, dataCandidates$Magnitude_1990)
# Confirmed dead and 2000 Magnitude
cor.test(dataCandidates$Confirmed_dead, dataCandidates$Magnitude_2000)
# Confirmed dead and Change in Magnitude
cor.test(dataCandidates$Confirmed_dead, dataCandidates$Magnitude_Change)

# log(Confirmed) dead and 1990 Magnitude
cor.test(dataCandidates$log_Confirmed_dead, dataCandidates$Magnitude_1990)
# log(Confirmed) and 2000 Magnitude
cor.test(dataCandidates$log_Confirmed_dead, dataCandidates$Magnitude_2000)
# log(Confirmed) and Change in Magnitude
cor.test(dataCandidates$log_Confirmed_dead, dataCandidates$Magnitude_Change)

##### OA6: Minimal Model Specifications #####
### NOTE ABOUT RESULTS FROM OA6: please note that the results from OA6 reflect one-tailed t-tests of directional hypotheses.
### Therefore, in order to obtain the correct p-values, (a) divide the p-values presented in the output by 2 and (b) remember
### that depending on which hypothesis is tested, only positive (Candidate Supply Hypothesis) or negative (Electoral Success 
### Hypothesis) coefficient estimates can be marked as statistically significant.

# Figure OA6.1: Bivariate Relationships between Dependent and Independent Variables
BasicModel1 <- lm(Candidates ~ log_Casualty, data=dataCandidates)
summary(BasicModel1) 
round(0.4695, 3)
# coefficient: 0.470
round(0.0518/2,3)
# p-value: 0.026
nobs(BasicModel1) # 87

BasicModel2 <- lm(Candidates ~ log_Confirmed_dead, data=dataCandidates)
summary(BasicModel2) 
round(0.5765, 3)
# coefficient: 0.577
round(0.026/2,3)
# p-value: 0.013
nobs(BasicModel2) # 87

BasicModel3 <- lm(Change_in_councilors ~ log_Casualty, data=dataCandidates)
summary(BasicModel3) 
round(-1.947, 3)
# coefficient: -1.947
round(0.0688/2,3)
# p-value: 0.034
nobs(BasicModel3) # 87

BasicModel4 <- lm(Change_in_councilors ~ log_Confirmed_dead, data=dataCandidates)
summary(BasicModel4) 
round(-2.2174, 3)
# coefficient: -2.217
round(0.0539/2,3)
# p-value: 0.027
nobs(BasicModel3) # 87

plot(dataCandidates$Candidates ~ dataCandidates$log_Casualty, pch=19, cex=0.5, xlab="log(Casualty)",
     ylab="Candidates", main="Bivariate Relationship between  \n Candidates and Casualty",
     ylim=c(25,40), xlim=c(-2.5, 2.5))
abline(BasicModel1, lty=2, lwd=2)
text(x = -1.25, y=34.5, substitute(paste(italic('coefficient'), " = 0.470")), cex=1, lwd=5, lwd=3)
text(x = -1.25, y=33.5, substitute(paste(italic('p-value'), " = 0.026")), cex=1, lwd=5, lwd=3)

plot(dataCandidates$Candidates ~ dataCandidates$log_Confirmed_dead, pch=19, cex=0.5, xlab="log(Confirmed dead)",
     ylab="Candidates", main="Bivariate Relationship between  \n Candidates and Confirmed Dead", 
     ylim=c(25,40), xlim=c(-2.5, 2.5 ))
abline(BasicModel2, lty=2, lwd=2)
text(x = -1.25, y=34.5, substitute(paste(italic('coefficient'), " = 0.577")), cex=1, lwd=5, lwd=3)
text(x = -1.25, y=33.5, substitute(paste(italic('p-value'), " = 0.013")), cex=1, lwd=5, lwd=2)

plot(dataCandidates$Change_in_councilors ~ dataCandidates$log_Casualty, pch=19, cex=0.5, xlab="log(Casualty)",
     ylab="Change in Councilors", main="Bivariate Relationship between \n Change in Councilors and Casualty", 
     ylim=c(-10, 40), xlim=c(-2.5, 2.5))
abline(BasicModel3, lty=2, lwd=2)
text(x = -1.5, y=11, substitute(paste(italic('coefficient'), " = -1.947")), cex=1, lwd=5, lwd=3)
text(x = -1.5, y=11-(10/3), substitute(paste(italic('p-value'), " = 0.034")), cex=1, lwd=5, lwd=2)

plot(dataCandidates$Change_in_councilors ~ dataCandidates$log_Confirmed_dead, pch=19, cex=0.5, xlab="log(Confirmed dead)",
     ylab="Change in Councilors", main="Bivariate Relationship between \n Change in Councilors and Confirmed Dead", 
     ylim=c(-10, 40), xlim=c(-2.5, 2.5))
abline(BasicModel4, lty=2, lwd=2)
text(x = -1.5, y=11, substitute(paste(italic('coefficient'), " = -2.217")), cex=1, lwd=5, lwd=3)
text(x = -1.5, y=11-(10/3), substitute(paste(italic('p-value'), " = 0.027")), cex=1, lwd=5, lwd=2)

# Table OA6.1: Wartime Violence and Female Candidates, Minimal Specifications
# Model 1: Candidates
summary(mod1 <- lm(Candidates ~ log_Casualty, data=dataCandidates))
nobs(mod1)

# Model 2: Candidates
summary(mod2 <- lm(Candidates ~ log_Confirmed_dead, data=dataCandidates))
nobs(mod2)

# Table OA6.2: Wartime Violence and Female Councilors, Minimal Specifications
# Model 1: Councilors
summary(mod1 <- lm(Change_in_councilors ~ log_Casualty, data=dataCandidates))
nobs(mod1)

# Model 2: Councilors
summary(mod2 <- lm(Change_in_councilors ~ log_Confirmed_dead, data=dataCandidates))
nobs(mod2)

##### OA8: Balance Tests #####
# Table OA8.1: Balance Tests for Casualty
# Model 1: Casualty
summary(mod1 <- lm(Casualty ~ Urban_share, data=dataCandidates))
nobs(mod1)
# Model 2: Casualty
summary(mod2 <- lm(Casualty ~ Income_per_capita, data=dataCandidates))
nobs(mod2)
# Model 3: Casualty
summary(mod3 <- lm(Casualty ~ Ethnic_polarization, data=dataCandidates))
nobs(mod3)
# Model 4: Casualty
summary(mod4 <- lm(Casualty ~ Share_of_women, data=dataCandidates))
nobs(mod4)
# Model 5: Casualty
summary(mod5 <- lm(Casualty ~ Y1990_councilors, data=dataCandidates))
nobs(mod5)
# Model 6: Casualty
summary(mod6 <- lm(Casualty ~ Urban_share + Income_per_capita + Ethnic_polarization +
                     Share_of_women + Y1990_councilors, data=dataCandidates))
nobs(mod6)

# Table OA8.2: Balance Tests for log(Casualty)
# Model 1: log(Casualty)
summary(mod1 <- lm(log_Casualty ~ Urban_share, data=dataCandidates))
nobs(mod1)
# Model 2: log(Casualty)
summary(mod2 <- lm(log_Casualty ~ Income_per_capita, data=dataCandidates))
nobs(mod2)
# Model 3: log(Casualty)
summary(mod3 <- lm(log_Casualty ~ Ethnic_polarization, data=dataCandidates))
nobs(mod3)
# Model 4: log(Casualty)
summary(mod4 <- lm(log_Casualty ~ Share_of_women, data=dataCandidates))
nobs(mod4)
# Model 5: log(Casualty)
summary(mod5 <- lm(log_Casualty ~ Y1990_councilors, data=dataCandidates))
nobs(mod5)
# Model 6: log(Casualty)
summary(mod6 <- lm(log_Casualty ~ Urban_share + Income_per_capita + Ethnic_polarization +
                     Share_of_women + Y1990_councilors, data=dataCandidates))
nobs(mod6)

# Table OA8.3: Balance Tests for Confirmed dead
# Model 1: Confirmed dead
summary(mod1 <- lm(Confirmed_dead ~ Urban_share, data=dataCandidates))
nobs(mod1)
# Model 2: Confirmed dead
summary(mod2 <- lm(Confirmed_dead ~ Income_per_capita, data=dataCandidates))
nobs(mod2)
# Model 3: Confirmed dead
summary(mod3 <- lm(Confirmed_dead ~ Ethnic_polarization, data=dataCandidates))
nobs(mod3)
# Model 4: Confirmed dead
summary(mod4 <- lm(Confirmed_dead ~ Share_of_women, data=dataCandidates))
nobs(mod4)
# Model 5: Confirmed dead
summary(mod5 <- lm(Confirmed_dead ~ Y1990_councilors, data=dataCandidates))
nobs(mod5)
# Model 6: Confirmed dead
summary(mod6 <- lm(Confirmed_dead ~ Urban_share + Income_per_capita + Ethnic_polarization +
                     Share_of_women + Y1990_councilors, data=dataCandidates))
nobs(mod6)

# Table OA8.4: Balance Tests for log(Confirmed dead)
# Model 1: log(Confirmed dead)
summary(mod1 <- lm(log_Confirmed_dead ~ Urban_share, data=dataCandidates))
nobs(mod1)
# Model 2: log(Confirmed dead)
summary(mod2 <- lm(log_Confirmed_dead ~ Income_per_capita, data=dataCandidates))
nobs(mod2)
# Model 3: log(Confirmed dead)
summary(mod3 <- lm(log_Confirmed_dead ~ Ethnic_polarization, data=dataCandidates))
nobs(mod3)
# Model 4: log(Confirmed dead)
summary(mod4 <- lm(log_Confirmed_dead ~ Share_of_women, data=dataCandidates))
nobs(mod4)
# Model 5: log(Confirmed dead)
summary(mod5 <- lm(log_Confirmed_dead ~ Y1990_councilors, data=dataCandidates))
nobs(mod5)
# Model 6: log(Confirmed dead)
summary(mod6 <- lm(log_Confirmed_dead ~ Urban_share + Income_per_capita + Ethnic_polarization +
                     Share_of_women + Y1990_councilors, data=dataCandidates))
nobs(mod6)

# Figure OA8.1: Variable Density Plots, Dichotomized at Median of log(Casualty)
data_treated <- dataCandidates[dataCandidates$log_Casualty > median(dataCandidates$log_Casualty),]
data_untreated <- dataCandidates[dataCandidates$log_Casualty <= median(dataCandidates$log_Casualty),]

wilcox.test(data_treated$Y1990_councilors, data_untreated$Y1990_councilors, exact=FALSE) 
plot(density(data_treated$Y1990_councilors), main='Councilors (in 1990)', xlab=substitute(paste(italic('p-value = 0.214'))))
par(new=TRUE)
plot(density(data_untreated$Y1990_councilors), lty=3, yaxt='n',xaxt='n', main='', xlab='')

wilcox.test(data_treated$Urban_share, data_untreated$Urban_share) 
plot(density(data_treated$Urban_share), main='Urban share', xlab=substitute(paste(italic('p-value = 0.460'))))
par(new=TRUE)
plot(density(data_untreated$Urban_share), lty=3, yaxt='n',xaxt='n', main='', xlab='')

wilcox.test(data_treated$Income_per_capita, data_untreated$Income_per_capita) 
plot(density(data_treated$Income_per_capita), main='Income per capita', xlab=substitute(paste(italic('p-value = 0.830'))))
par(new=TRUE)
plot(density(data_untreated$Income_per_capita), lty=3, yaxt='n',xaxt='n', main='', xlab='')

wilcox.test(data_treated$Ethnic_polarization, data_untreated$Ethnic_polarization) # 0.000
plot(density(data_treated$Ethnic_polarization), main='Ethnic polarization', xlab=substitute(paste(italic('p-value = 0.000'))))
par(new=TRUE)
plot(density(data_untreated$Ethnic_polarization), lty=3, yaxt='n',xaxt='n', main='', xlab='')

wilcox.test(data_treated$Share_of_women, data_untreated$Share_of_women) 
plot(density(data_treated$Share_of_women), main='Share of women', xlab=substitute(paste(italic('p-value = 0.676'))))
par(new=TRUE)
plot(density(data_untreated$Share_of_women), lty=3, yaxt='n',xaxt='n', main='', xlab='')

plot.new()
legend('center', lty=c(1,3), legend = c('Above median', 'At or below median '), bty='n')

# Figure OA8.2: Variable Density Plots, Dichotomized at Median of log(Confirmed dead)
data_treated <- dataCandidates[dataCandidates$log_Confirmed_dead > median(dataCandidates$log_Confirmed_dead),]
data_untreated <- dataCandidates[dataCandidates$log_Confirmed_dead <= median(dataCandidates$log_Confirmed_dead),]

wilcox.test(data_treated$Y1990_councilors, data_untreated$Y1990_councilors, exact=FALSE) 
plot(density(data_treated$Y1990_councilors), main='Councilors (in 1990)', xlab=substitute(paste(italic('p-value = 0.469'))))
par(new=TRUE)
plot(density(data_untreated$Y1990_councilors), lty=3, yaxt='n',xaxt='n', main='', xlab='')

wilcox.test(data_treated$Urban_share, data_untreated$Urban_share)
plot(density(data_treated$Urban_share), main='Urban share', xlab=substitute(paste(italic('p-value = 0.424'))))
par(new=TRUE)
plot(density(data_untreated$Urban_share), lty=3, yaxt='n',xaxt='n', main='', xlab='')

wilcox.test(data_treated$Income_per_capita, data_untreated$Income_per_capita) 
plot(density(data_treated$Income_per_capita), main='Income per capita', xlab=substitute(paste(italic('p-value = 0.950'))))
par(new=TRUE)
plot(density(data_untreated$Income_per_capita), lty=3, yaxt='n',xaxt='n', main='', xlab='')

wilcox.test(data_treated$Ethnic_polarization, data_untreated$Ethnic_polarization) 
plot(density(data_treated$Ethnic_polarization), main='Ethnic polarization', xlab=substitute(paste(italic('p-value = 0.000'))))
par(new=TRUE)
plot(density(data_untreated$Ethnic_polarization), lty=3, yaxt='n',xaxt='n', main='', xlab='')

wilcox.test(data_treated$Share_of_women, data_untreated$Share_of_women) 
plot(density(data_treated$Share_of_women), main='Share of women', xlab=substitute(paste(italic('p-value = 0.790'))))
par(new=TRUE)
plot(density(data_untreated$Share_of_women), lty=3, yaxt='n',xaxt='n', main='', xlab='')

plot.new()
legend('center', lty=c(1,3), legend = c('Above median', 'At or below median '), bty='n')

##### OA9: Supplementary Model Results #####
### NOTE ABOUT RESULTS FROM OA9: please note that the results from OA9 in some instances reflect two-tailed t-tests and in other
### instances one-tailed t-tests of directional hypotheses. Please consult the online appendix in order to see when two-tailed
### and one-tailed t-tests are used (it is noted at the bottom of every table). In cases of one-tailed t-tests, in order to obtain 
### the correct p-values, (a) divide the p-values presented in the output by 2 and (b) remember that depending on which hypothesis 
### is tested, only positive or negative coefficient estimates can be marked as statistically significant.

# Table OA9.1: Wartime Violence and Female Candidates, Population Density Control
# Model 1: Candidates
summary(lm(Candidates ~ log_Casualty + Population_density + Income_per_capita +
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 2: Candidates
summary(lm(Candidates ~ log_Confirmed_dead + Population_density + Income_per_capita +
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))

# Table OA9.2: Wartime Violence and Female Candidates, Population Size Control
# Model 1: Candidates
summary(lm(Candidates ~ log_Casualty + log(Population_size) + Income_per_capita +
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 2: Candidates
summary(lm(Candidates ~ log_Confirmed_dead + log(Population_size) + Income_per_capita +
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))

# Table OA9.3: Wartime Violence and Female Candidates, Ethnic Fractionalization Control
# Model 1: Candidates
summary(lm(Candidates ~ log_Casualty + Urban_share + Income_per_capita +
             Ethnic_fractionalization + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 2: Candidates
summary(lm(Candidates ~ log_Confirmed_dead + Urban_share + Income_per_capita +
             Ethnic_fractionalization + Share_of_women + Y1990_councilors, data=dataCandidates))

# Table OA9.4: Wartime Violence and Female Candidates, Ethnic Polarization and Fractionalization Controls
# Model 1: Candidates
summary(lm(Candidates ~ log_Casualty + Urban_share + Income_per_capita + Ethnic_polarization +
             Ethnic_fractionalization + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 2: Candidates
summary(lm(Candidates ~ log_Confirmed_dead + Urban_share + Income_per_capita + Ethnic_polarization +
             Ethnic_fractionalization + Share_of_women + Y1990_councilors, data=dataCandidates))

# Table OA9.5: Wartime Violence and Female Candidates, Ethnic Majority/Plurality Dummy Controls
# Model 1: Candidates
summary(lm(Candidates ~ log_Casualty + Urban_share + Income_per_capita  +
             Majority_Plurality + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 2: Candidates
summary(lm(Candidates ~ log_Confirmed_dead + Urban_share + Income_per_capita +
             Majority_Plurality + Share_of_women + Y1990_councilors, data=dataCandidates))

dataCandidates$Majority <- relevel(dataCandidates$Majority, ref="No majority group")
# Table OA9.6: Wartime Violence and Female Candidates, Ethnic Majority Dummy Controls
# Model 1: Candidates
summary(lm(Candidates ~ log_Casualty + Urban_share + Income_per_capita  +
             Majority + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 2: Candidates
summary(lm(Candidates ~ log_Confirmed_dead + Urban_share + Income_per_capita +
             Majority + Share_of_women + Y1990_councilors, data=dataCandidates))

# Table OA9.7: Wartime Violence and Female Candidates, Economic and Social Gender Inequality Controls
# Model 1: Candidates
summary(lm(Candidates ~ log_Casualty + Urban_share + Income_per_capita + Ethnic_polarization +
             Share_of_women + Y1990_councilors + Employment_share + Illiteracy_share, data=dataCandidates))
# Model 2: Candidates
summary(lm(Candidates ~ log_Confirmed_dead + Urban_share + Income_per_capita + Ethnic_polarization +
             Share_of_women + Y1990_councilors + Employment_share + Illiteracy_share, data=dataCandidates))

dataCandidates$Candidates_proportion <- (dataCandidates$Candidates)/100
# Table OA9.8: Wartime Violence and Female Candidates, Beta Regression
# Model 1: Candidates proportion
summary(betareg(Candidates_proportion ~ log_Casualty + Urban_share + Income_per_capita + 
                  Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 2: Candidates proportion
summary(betareg(Candidates_proportion ~ log_Confirmed_dead + Urban_share + Income_per_capita + 
                  Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))

dataCandidates$Logit_transformed_candidate <- logit(dataCandidates$Candidates_proportion)
# Table OA9.9: Wartime Violence and Female Candidates, Logit Transformed Candidates
summary(lm(Logit_transformed_candidate ~ log_Casualty + Urban_share + Income_per_capita + 
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))
summary(lm(Logit_transformed_candidate ~ log_Confirmed_dead + Urban_share + Income_per_capita + 
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))

# Table OA9.10: Wartime Violence and Female Councilors, Two-Limit Tobit Regression
# Model 1: Councilors
summary(censReg(Councilors ~ log_Casualty:PrePost + PrePost + Municipality - 1, left=0, right=100, data=dataCouncilors))
# Model 2: Councilors
summary(censReg(Councilors ~ log_Confirmed_dead:PrePost + PrePost + Municipality - 1, left=0, right=100, data=dataCouncilors))

# Table OA9.11: Wartime Violence and the Number of Male and Female Candidates, Negative Binomial Models
# Model 1: Female candidates
summary(glm.nb(Female_candidates ~ log_Casualty + Urban_share + Income_per_capita + 
                 Ethnic_polarization + Share_of_women + Y1990_councilors +
                 log(Population_size), data=dataCandidates))
# Model 2: Male candidates
summary(glm.nb(Male_candidates ~ log_Casualty + Urban_share + Income_per_capita + 
                 Ethnic_polarization + Share_of_women + Y1990_councilors +
                 log(Population_size), data=dataCandidates))
# Model 3: Female candidates
summary(glm.nb(Female_candidates ~ log_Confirmed_dead + Urban_share + Income_per_capita + 
                 Ethnic_polarization + Share_of_women + Y1990_councilors +
                 log(Population_size), data=dataCandidates))
# Model 4: Male candidates
summary(glm.nb(Male_candidates ~ log_Confirmed_dead + Urban_share + Income_per_capita + 
                 Ethnic_polarization + Share_of_women + Y1990_councilors +
                 log(Population_size), data=dataCandidates))

dataMilitary <- dataCouncilors[complete.cases(dataCouncilors$Military_share),]
# Table OA9.12: Wartime Violence and Female Councilors, Military Measures
# Model 1: Councilors
m1 <-  lm(Councilors ~ Military_share:PrePost + PrePost + Municipality - 1, data=dataMilitary)
vcovCL1 <- cluster.vcov(m1, ~dataMilitary$Municipality+dataMilitary$PrePost)
coeftest(m1, vcovCL1) 
nobs(m1) 
summary(m1)$r.squared 
# Model 2: Councilors
m2 <-  lm(Councilors ~ log_Military_share:PrePost + PrePost + Municipality - 1, data=dataMilitary)
vcovCL1 <- cluster.vcov(m2, ~dataMilitary$Municipality+dataMilitary$PrePost)
coeftest(m2, vcovCL1)
nobs(m2) 
summary(m2)$r.squared 
names(dataMilitary)
# Model 3: Councilors
m3 <-  lm(Councilors ~ Casualty_difference:PrePost + PrePost + Municipality - 1, data=dataMilitary)
vcovCL1 <- cluster.vcov(m3, ~dataMilitary$Municipality+dataMilitary$PrePost)
coeftest(m3, vcovCL1)
nobs(m3) 
summary(m3)$r.squared

dataCandidates$Male_vote_fractionalization <- (dataCandidates$Male_vote_fractionalization)*100
dataCandidates$Female_vote_fractionalization <- (dataCandidates$Female_vote_fractionalization)*100
# Table OA9.13: Wartime Violence and Male Vote Fractionalization
# Model 1: Male vote fractionalization
summary(lm(Male_vote_fractionalization ~ log_Casualty + Urban_share + Income_per_capita + Ethnic_polarization +
             Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 2: Male vote fractionalization
summary(lm(Male_vote_fractionalization ~ log_Confirmed_dead + Urban_share + Income_per_capita + Ethnic_polarization +
             Share_of_women + Y1990_councilors, data=dataCandidates))

# Table OA9.14: Wartime Violence and Female Vote Fractionalization
# Model 1: Female vote fractionalization
summary(lm(Female_vote_fractionalization ~ log_Casualty + Urban_share + Income_per_capita + Ethnic_polarization +
             Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 2: Female vote fractionalization
summary(lm(Female_vote_fractionalization ~ log_Confirmed_dead + Urban_share + Income_per_capita + Ethnic_polarization +
             Share_of_women + Y1990_councilors, data=dataCandidates))

# Table OA9.15: Wartime Violence and Female Candidates, Mediating Effect of Gender Balance
# Mode 1: Candidates
first <-  lm(Candidates ~ log_Casualty + Urban_share + Income_per_capita + Ethnic_polarization +
               Share_of_women + Y1990_councilors + Share_of_women_2013, data=dataCandidates)

direct <- lm(I(Candidates - coef(first)['Share_of_women_2013']*Share_of_women_2013) ~ log_Casualty + 
               Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + 
               Y1990_councilors, data=dataCandidates)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=7)
for(b in 1:boots){
  d.star <- dataCandidates[sample(1:nrow(dataCandidates), replace=TRUE),]
  first.boots <- lm(Candidates ~ log_Casualty + Urban_share + Income_per_capita + Ethnic_polarization +
                      Share_of_women + Y1990_councilors + Share_of_women_2013, d.star)
  direct.boots <- lm(I(Candidates - coef(first.boots)['Share_of_women_2013']*Share_of_women_2013) ~ log_Casualty + 
                       Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + Y1990_councilors, d.star)
  fl.boots[b,] <- coef(direct.boots)
}
# 80 degrees of freedom
summary(direct)$r.squared
coef(direct) # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.95), df=80) # t-value: 1.664125

# Mode 2: Candidates
first <-  lm(Candidates ~ log_Confirmed_dead + Urban_share + Income_per_capita + Ethnic_polarization +
               Share_of_women + Y1990_councilors + Share_of_women_2013, data=dataCandidates)
direct <- lm(I(Candidates - coef(first)['Share_of_women_2013']*Share_of_women_2013) ~ log_Confirmed_dead + 
               Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + 
               Y1990_councilors, data=dataCandidates)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=7)
for(b in 1:boots){
  d.star <- dataCandidates[sample(1:nrow(dataCandidates), replace=TRUE),]
  first.boots <- lm(Candidates ~ log_Confirmed_dead + Urban_share + Income_per_capita + Ethnic_polarization +
                      Share_of_women + Y1990_councilors + Share_of_women_2013, d.star)
  direct.boots <- lm(I(Candidates - coef(first.boots)['Share_of_women_2013']*Share_of_women_2013) ~ log_Confirmed_dead + 
                       Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + Y1990_councilors, d.star)
  fl.boots[b,] <- coef(direct.boots)
}
# 80 degrees of freedom
summary(direct)$r.squared
coef(direct) # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.95), df=80) # t-value: 1.664125

# Table OA9.16: Wartime Violence and Female Councilors, Mediating Effect of Gender Balance
# Model 1: Councilors
first <-  lm(Councilors ~ log_Casualty:PrePost + PrePost + Municipality - 1 + Share_of_women, data=dataCouncilors)
direct <- lm(I(Councilors - coef(first)['Share_of_women']*Share_of_women) ~ log_Casualty:PrePost + PrePost + Municipality - 1, data=dataCouncilors)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=2)
for(b in 1:boots){
  d.star <- dataCouncilors[sample(1:nrow(dataCouncilors), replace=TRUE),]
  first.boots <- lm(Councilors ~ log_Casualty:PrePost + PrePost + Share_of_women + Municipality - 1, d.star)
  direct.boots <- lm(I(Councilors - coef(first.boots)['Share_of_women']*Share_of_women) ~ log_Casualty:PrePost + PrePost + Municipality - 1, d.star)
  fl.boots[b,] <- coef(direct.boots)[c(1, (length(coef(direct.boots))):length(coef(direct.boots)))]
}
# 85 degrees of freedom
summary(direct)$r.squared
coef(direct)[c(1,89)] # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.05), df=85) # t-value: -1.662978

# Model 2: Councilors
first <-  lm(Councilors ~ log_Confirmed_dead:PrePost + PrePost + Municipality - 1 + Share_of_women, data=dataCouncilors)
direct <- lm(I(Councilors - coef(first)['Share_of_women']*Share_of_women) ~ log_Confirmed_dead:PrePost + PrePost + Municipality - 1, data=dataCouncilors)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=2)
for(b in 1:boots){
  d.star <- dataCouncilors[sample(1:nrow(dataCouncilors), replace=TRUE),]
  first.boots <- lm(Councilors ~ log_Confirmed_dead:PrePost + PrePost + Share_of_women + Municipality - 1, d.star)
  direct.boots <- lm(I(Councilors - coef(first.boots)['Share_of_women']*Share_of_women) ~ log_Confirmed_dead:PrePost + PrePost + Municipality - 1, d.star)
  fl.boots[b,] <- coef(direct.boots)[c(1, (length(coef(direct.boots))):length(coef(direct.boots)))]
}
# 85 degrees of freedom
summary(direct)$r.squared
coef(direct)[c(1,89)] # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.05), df=85) # t-value: -1.662978

# Table OA9.17: Wartime Violence and Women's List Placement, Relative Placement Outcome
# Model 1: Relative placement
summary(lm(Relative_placement ~ log_Casualty + Urban_share + Income_per_capita +
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 2: Relative placement
summary(lm(Relative_placement ~ log_Confirmed_dead + Urban_share + Income_per_capita +
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 3: Relative placement
first <-  lm(Relative_placement ~ log_Casualty + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors + Candidates, data=dataCandidates)
direct <- lm(I(Relative_placement - coef(first)['Candidates']*Candidates) ~ log_Casualty + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=7)
for(b in 1:boots){
  d.star <- dataCandidates[sample(1:nrow(dataCandidates), replace=TRUE),]
  first.boots <- lm(Relative_placement ~ log_Casualty + Urban_share + Income_per_capita + Ethnic_polarization +
                      Share_of_women + Y1990_councilors + Candidates, d.star)
  direct.boots <- lm(I(Relative_placement - coef(first.boots)['Candidates']*Candidates) ~ log_Casualty + 
                       Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + Y1990_councilors, d.star)
  fl.boots[b,] <- coef(direct.boots)
}
# 80 degrees of freedom
summary(direct)$r.squared
coef(direct) # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.025, 0.975), df=80) # t-value: -1.990063 to 1.990063

# Model 4: Relative placement
first <-  lm(Relative_placement ~ log_Confirmed_dead + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors + Candidates, data=dataCandidates)
direct <- lm(I(Relative_placement - coef(first)['Candidates']*Candidates) ~ log_Confirmed_dead + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=7)
for(b in 1:boots){
  d.star <- dataCandidates[sample(1:nrow(dataCandidates), replace=TRUE),]
  first.boots <- lm(Relative_placement ~ log_Confirmed_dead + Urban_share + Income_per_capita + Ethnic_polarization +
                      Share_of_women + Y1990_councilors + Candidates, d.star)
  direct.boots <- lm(I(Relative_placement - coef(first.boots)['Candidates']*Candidates) ~ log_Confirmed_dead + 
                       Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + Y1990_councilors, d.star)
  fl.boots[b,] <- coef(direct.boots)
}
# 80 degrees of freedom
summary(direct)$r.squared
coef(direct) # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.025, 0.975), df=80) # t-value: -1.990063 to 1.990063

# Table OA9.18: Wartime Violence and Women's List Placement, Top 10 Outcome
# Model 1: Top 10
summary(lm(Top_10 ~ log_Casualty + Urban_share + Income_per_capita +
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 2: Top 10
summary(lm(Top_10 ~ log_Confirmed_dead + Urban_share + Income_per_capita +
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 3: Top 10
first <-  lm(Top_10 ~ log_Casualty + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors + Candidates, data=dataCandidates)
direct <- lm(I(Top_10 - coef(first)['Candidates']*Candidates) ~ log_Casualty + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=7)
for(b in 1:boots){
  d.star <- dataCandidates[sample(1:nrow(dataCandidates), replace=TRUE),]
  first.boots <- lm(Top_10 ~ log_Casualty + Urban_share + Income_per_capita + Ethnic_polarization +
                      Share_of_women + Y1990_councilors + Candidates, d.star)
  direct.boots <- lm(I(Top_10 - coef(first.boots)['Candidates']*Candidates) ~ log_Casualty + 
                       Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + Y1990_councilors, d.star)
  fl.boots[b,] <- coef(direct.boots)
}
# 80 degrees of freedom
summary(direct)$r.squared
coef(direct) # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.025, 0.975), df=80) # t-value: -1.990063 to 1.990063

# Model 4: Top 10
first <-  lm(Top_10 ~ log_Confirmed_dead + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors + Candidates, data=dataCandidates)
direct <- lm(I(Top_10 - coef(first)['Candidates']*Candidates) ~ log_Confirmed_dead + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=7)
for(b in 1:boots){
  d.star <- dataCandidates[sample(1:nrow(dataCandidates), replace=TRUE),]
  first.boots <- lm(Top_10 ~ log_Confirmed_dead + Urban_share + Income_per_capita + Ethnic_polarization +
                      Share_of_women + Y1990_councilors + Candidates, d.star)
  direct.boots <- lm(I(Top_10 - coef(first.boots)['Candidates']*Candidates) ~ log_Confirmed_dead + 
                       Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + Y1990_councilors, d.star)
  fl.boots[b,] <- coef(direct.boots)
}
# 80 degrees of freedom
summary(direct)$r.squared
coef(direct) # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.025, 0.975), df=80) # t-value: -1.990063 to 1.990063

# Table OA9.19: Wartime Violence and Women's List Placement, Top 5 Outcome
# Model 1: Top 5
summary(lm(Top_5 ~ log_Casualty + Urban_share + Income_per_capita +
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 2: Top 5
summary(lm(Top_5 ~ log_Confirmed_dead + Urban_share + Income_per_capita +
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 3: Top 5
first <-  lm(Top_5 ~ log_Casualty + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors + Candidates, data=dataCandidates)
direct <- lm(I(Top_5 - coef(first)['Candidates']*Candidates) ~ log_Casualty + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=7)
for(b in 1:boots){
  d.star <- dataCandidates[sample(1:nrow(dataCandidates), replace=TRUE),]
  first.boots <- lm(Top_5 ~ log_Casualty + Urban_share + Income_per_capita + Ethnic_polarization +
                      Share_of_women + Y1990_councilors + Candidates, d.star)
  direct.boots <- lm(I(Top_5 - coef(first.boots)['Candidates']*Candidates) ~ log_Casualty + 
                       Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + Y1990_councilors, d.star)
  fl.boots[b,] <- coef(direct.boots)
}
# 80 degrees of freedom
summary(direct)$r.squared
coef(direct) # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.025, 0.975), df=80) # t-value: -1.990063 to 1.990063

# Model 4: Top 5
first <-  lm(Top_5 ~ log_Confirmed_dead + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors + Candidates, data=dataCandidates)
direct <- lm(I(Top_5 - coef(first)['Candidates']*Candidates) ~ log_Confirmed_dead + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=7)
for(b in 1:boots){
  d.star <- dataCandidates[sample(1:nrow(dataCandidates), replace=TRUE),]
  first.boots <- lm(Top_5 ~ log_Confirmed_dead + Urban_share + Income_per_capita + Ethnic_polarization +
                      Share_of_women + Y1990_councilors + Candidates, d.star)
  direct.boots <- lm(I(Top_5 - coef(first.boots)['Candidates']*Candidates) ~ log_Confirmed_dead + 
                       Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + Y1990_councilors, d.star)
  fl.boots[b,] <- coef(direct.boots)
}
# 80 degrees of freedom
summary(direct)$r.squared
coef(direct) # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.025, 0.975), df=80) # t-value: -1.990063 to 1.990063

# Table OA9.20: Wartime Violence and Women's List Placement, Top 4 Outcome
# Model 1: Top 4
summary(lm(Top_4 ~ log_Casualty + Urban_share + Income_per_capita +
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 2: Top 4
summary(lm(Top_4 ~ log_Confirmed_dead + Urban_share + Income_per_capita +
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 3: Top 4
first <-  lm(Top_4 ~ log_Casualty + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors + Candidates, data=dataCandidates)
direct <- lm(I(Top_4 - coef(first)['Candidates']*Candidates) ~ log_Casualty + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=7)
for(b in 1:boots){
  d.star <- dataCandidates[sample(1:nrow(dataCandidates), replace=TRUE),]
  first.boots <- lm(Top_4 ~ log_Casualty + Urban_share + Income_per_capita + Ethnic_polarization +
                      Share_of_women + Y1990_councilors + Candidates, d.star)
  direct.boots <- lm(I(Top_4 - coef(first.boots)['Candidates']*Candidates) ~ log_Casualty + 
                       Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + Y1990_councilors, d.star)
  fl.boots[b,] <- coef(direct.boots)
}
# 80 degrees of freedom
summary(direct)$r.squared
coef(direct) # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.025, 0.975), df=80) # t-value: -1.990063 to 1.990063

# Model 4: Top 4
first <-  lm(Top_4 ~ log_Confirmed_dead + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors + Candidates, data=dataCandidates)
direct <- lm(I(Top_4 - coef(first)['Candidates']*Candidates) ~ log_Confirmed_dead + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=7)
for(b in 1:boots){
  d.star <- dataCandidates[sample(1:nrow(dataCandidates), replace=TRUE),]
  first.boots <- lm(Top_4 ~ log_Confirmed_dead + Urban_share + Income_per_capita + Ethnic_polarization +
                      Share_of_women + Y1990_councilors + Candidates, d.star)
  direct.boots <- lm(I(Top_4 - coef(first.boots)['Candidates']*Candidates) ~ log_Confirmed_dead + 
                       Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + Y1990_councilors, d.star)
  fl.boots[b,] <- coef(direct.boots)
}
# 80 degrees of freedom
summary(direct)$r.squared
coef(direct) # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.025, 0.975), df=80) # t-value: -1.990063 to 1.990063

# Table OA9.21: Wartime Violence and Women's List Placement, Top 3 Outcome
# Model 1: Top 3
summary(lm(Top_3 ~ log_Casualty + Urban_share + Income_per_capita +
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 2: Top 3
summary(lm(Top_3 ~ log_Confirmed_dead + Urban_share + Income_per_capita +
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 3: Top 3
first <-  lm(Top_3 ~ log_Casualty + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors + Candidates, data=dataCandidates)
direct <- lm(I(Top_3 - coef(first)['Candidates']*Candidates) ~ log_Casualty + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=7)
for(b in 1:boots){
  d.star <- dataCandidates[sample(1:nrow(dataCandidates), replace=TRUE),]
  first.boots <- lm(Top_3 ~ log_Casualty + Urban_share + Income_per_capita + Ethnic_polarization +
                      Share_of_women + Y1990_councilors + Candidates, d.star)
  direct.boots <- lm(I(Top_3 - coef(first.boots)['Candidates']*Candidates) ~ log_Casualty + 
                       Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + Y1990_councilors, d.star)
  fl.boots[b,] <- coef(direct.boots)
}
# 80 degrees of freedom
summary(direct)$r.squared
coef(direct) # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.025, 0.975), df=80) # t-value: -1.990063 to 1.990063

# Model 4: Top 3
first <-  lm(Top_3 ~ log_Confirmed_dead + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors + Candidates, data=dataCandidates)
direct <- lm(I(Top_3 - coef(first)['Candidates']*Candidates) ~ log_Confirmed_dead + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=7)
for(b in 1:boots){
  d.star <- dataCandidates[sample(1:nrow(dataCandidates), replace=TRUE),]
  first.boots <- lm(Top_3 ~ log_Confirmed_dead + Urban_share + Income_per_capita + Ethnic_polarization +
                      Share_of_women + Y1990_councilors + Candidates, d.star)
  direct.boots <- lm(I(Top_3 - coef(first.boots)['Candidates']*Candidates) ~ log_Confirmed_dead + 
                       Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + Y1990_councilors, d.star)
  fl.boots[b,] <- coef(direct.boots)
}
# 80 degrees of freedom
summary(direct)$r.squared
coef(direct) # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.025, 0.975), df=80) # t-value: -1.990063 to 1.990063

# Table OA9.22: Wartime Violence and Women's List Placement, Top 2 Outcome
# Model 1: Top 2
summary(lm(Top_2 ~ log_Casualty + Urban_share + Income_per_capita +
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 2: Top 2
summary(lm(Top_2 ~ log_Confirmed_dead + Urban_share + Income_per_capita +
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 3: Top 2
first <-  lm(Top_2 ~ log_Casualty + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors + Candidates, data=dataCandidates)
direct <- lm(I(Top_2 - coef(first)['Candidates']*Candidates) ~ log_Casualty + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=7)
for(b in 1:boots){
  d.star <- dataCandidates[sample(1:nrow(dataCandidates), replace=TRUE),]
  first.boots <- lm(Top_2 ~ log_Casualty + Urban_share + Income_per_capita + Ethnic_polarization +
                      Share_of_women + Y1990_councilors + Candidates, d.star)
  direct.boots <- lm(I(Top_2 - coef(first.boots)['Candidates']*Candidates) ~ log_Casualty + 
                       Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + Y1990_councilors, d.star)
  fl.boots[b,] <- coef(direct.boots)
}
# 80 degrees of freedom
summary(direct)$r.squared
coef(direct) # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.025, 0.975), df=80) # t-value: -1.990063 to 1.990063

# Model 4: Top 2
first <-  lm(Top_2 ~ log_Confirmed_dead + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors + Candidates, data=dataCandidates)
direct <- lm(I(Top_2 - coef(first)['Candidates']*Candidates) ~ log_Confirmed_dead + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=7)
for(b in 1:boots){
  d.star <- dataCandidates[sample(1:nrow(dataCandidates), replace=TRUE),]
  first.boots <- lm(Top_2 ~ log_Confirmed_dead + Urban_share + Income_per_capita + Ethnic_polarization +
                      Share_of_women + Y1990_councilors + Candidates, d.star)
  direct.boots <- lm(I(Top_2 - coef(first.boots)['Candidates']*Candidates) ~ log_Confirmed_dead + 
                       Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + Y1990_councilors, d.star)
  fl.boots[b,] <- coef(direct.boots)
}
# 80 degrees of freedom
summary(direct)$r.squared
coef(direct) # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.025, 0.975), df=80) # t-value: -1.990063 to 1.990063

# Table OA9.23: Wartime Violence and Women's List Placement, Top 1 Outcome
# Model 1: Top 1
summary(lm(Top_1 ~ log_Casualty + Urban_share + Income_per_capita +
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 2: Top 1
summary(lm(Top_1 ~ log_Confirmed_dead + Urban_share + Income_per_capita +
             Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates))
# Model 3: Top 1
first <-  lm(Top_1 ~ log_Casualty + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors + Candidates, data=dataCandidates)
direct <- lm(I(Top_1 - coef(first)['Candidates']*Candidates) ~ log_Casualty + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=7)
for(b in 1:boots){
  d.star <- dataCandidates[sample(1:nrow(dataCandidates), replace=TRUE),]
  first.boots <- lm(Top_1 ~ log_Casualty + Urban_share + Income_per_capita + Ethnic_polarization +
                      Share_of_women + Y1990_councilors + Candidates, d.star)
  direct.boots <- lm(I(Top_1 - coef(first.boots)['Candidates']*Candidates) ~ log_Casualty + 
                       Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + Y1990_councilors, d.star)
  fl.boots[b,] <- coef(direct.boots)
}
# 80 degrees of freedom
summary(direct)$r.squared
coef(direct) # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.025, 0.975), df=80) # t-value: -1.990063 to 1.990063

# Model 4: Top 1
first <-  lm(Top_1 ~ log_Confirmed_dead + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors + Candidates, data=dataCandidates)
direct <- lm(I(Top_1 - coef(first)['Candidates']*Candidates) ~ log_Confirmed_dead + Urban_share + Income_per_capita +
               Ethnic_polarization + Share_of_women + Y1990_councilors, data=dataCandidates)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=7)
for(b in 1:boots){
  d.star <- dataCandidates[sample(1:nrow(dataCandidates), replace=TRUE),]
  first.boots <- lm(Top_1 ~ log_Confirmed_dead + Urban_share + Income_per_capita + Ethnic_polarization +
                      Share_of_women + Y1990_councilors + Candidates, d.star)
  direct.boots <- lm(I(Top_1 - coef(first.boots)['Candidates']*Candidates) ~ log_Confirmed_dead + 
                       Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + Y1990_councilors, d.star)
  fl.boots[b,] <- coef(direct.boots)
}
# 80 degrees of freedom
summary(direct)$r.squared
coef(direct) # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.025, 0.975), df=80) # t-value: -1.990063 to 1.990063

# Table OA9.24: Wartime Violence and Female Candidates, Mediating Effect of Average List Length
# Model 1: Candidates
first <-  lm(Candidates ~ log_Casualty + Urban_share + Income_per_capita + Ethnic_polarization +
               Share_of_women + Y1990_councilors + Average_list_length, data=dataCandidates)
direct <- lm(I(Candidates - coef(first)['Average_list_length']*Average_list_length) ~ log_Casualty + 
               Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + 
               Y1990_councilors, data=dataCandidates)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=7)
for(b in 1:boots){
  d.star <- dataCandidates[sample(1:nrow(dataCandidates), replace=TRUE),]
  first.boots <- lm(Candidates ~ log_Casualty + Urban_share + Income_per_capita + Ethnic_polarization +
                      Share_of_women + Y1990_councilors + Average_list_length, d.star)
  direct.boots <- lm(I(Candidates - coef(first.boots)['Average_list_length']*Average_list_length) ~ log_Casualty + 
                       Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + Y1990_councilors, d.star)
  fl.boots[b,] <- coef(direct.boots)
}
# 80 degrees of freedom
summary(direct)$r.squared
coef(direct) # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.95), df=80) # t-value: 1.664125

# Model 2: Candidates
first <-  lm(Candidates ~ log_Confirmed_dead + Urban_share + Income_per_capita + Ethnic_polarization +
               Share_of_women + Y1990_councilors + Average_list_length, data=dataCandidates)
direct <- lm(I(Candidates - coef(first)['Average_list_length']*Average_list_length) ~ log_Confirmed_dead + 
               Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + 
               Y1990_councilors, data=dataCandidates)
set.seed(1)
boots <- 1000
fl.boots <- matrix(NA, nrow=boots, ncol=7)
for(b in 1:boots){
  d.star <- dataCandidates[sample(1:nrow(dataCandidates), replace=TRUE),]
  first.boots <- lm(Candidates ~ log_Confirmed_dead + Urban_share + Income_per_capita + Ethnic_polarization +
                      Share_of_women + Y1990_councilors + Average_list_length, d.star)
  direct.boots <- lm(I(Candidates - coef(first.boots)['Average_list_length']*Average_list_length) ~ log_Confirmed_dead + 
                       Urban_share + Income_per_capita + Ethnic_polarization + Share_of_women + Y1990_councilors, d.star)
  fl.boots[b,] <- coef(direct.boots)
}
# 80 degrees of freedom
summary(direct)$r.squared
coef(direct) # coefficients
(SEs <- apply(fl.boots,2,sd)) # standard errors
qt(c(0.95), df=80) # t-value: 1.664125
